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We address several questions on the behavior of a numerical model recently introduced to study 
seismic phenomena, that includes relaxation in the plates as a key ingredient. We make an analysis of 
the scaling of the largest events with system size, and show that when parameters are appropriately 
interpreted, the typical size of the largest events scale as the system size, without the necessity 
to tune any parameter. Secondly, we show that the temporal activity in the model is inherently 
non-stationary, and obtain from here justification and support for the concept of a "seismic cycle" 
in the temporal evolution of seismic activity. Finally, we ask for the reasons that make the model 
display a realistic value of the decaying exponent h in the Gutenberg-Richter law for the avalanche 
size distribution. We explain why relaxation induces a systematic increase of h from its value 
b ~ 0.4 observed in the absence of relaxation. However, we have not been able to justify the actual 
robustness of the model in displaying a consistent b value around the experimentally observed value 
6~ 1. 



I. INTRODUCTION 



Modeling of seismic phenomena as a statistical me- 
chanic process has a long history that goes back almost 
to the very beginning of plate tectonics. There are a lot 
of variables affecting seismic activity, and many of them 
cannot be taken into account in a simplified theoreti- 
cal description. However, if a model gives a sequence of 
events strongly resembling real ones (in its temporal, spa- 
tial, and magnitude distribution), it is tempting to con- 
sider that the model captures those key ingredients play- 
ing an important role in the production of seismic events. 
This is what happens with the model proposed by Bur- 
ridge and Knopoffi^. Constructed in the form of a large 
collection of repetitive and simple elements (rigid blocks 
interacting through elastic springs and sliding on a rigid 
surface), it displays non-trivial features, most remarkably 
the existence of avalanches (earthquakes) with a broad 
distribution of sizes. Based on the work of Burridge and 
Knopoff, in 1991 Olami, Feder and Christensen^ (OFC) 
introduced a model that has become one of the paradigms 
of simulation of seismic activity. The model is presented 
in the form of cellular automata, and the dynamics is 
defined as a list of rules. A broad distribution of event 
sizes is obtained with the OFC model. Despite the at- 
tention it received, the OFC model displays a series of 
unrealistic features when trying to model seismic phe- 
nomena. For instance, the prominent spatial and tem- 
poral correlation of events that is observed in the field, 
mostly the existence of aftershocks, is not reproduced^. 
Also, obtaining a realistic value of the b exponent in the 
Gutenberg-Richter (GR) law requires the tuning of an 
internal parameter of the model. 

In the last years, we have been working in introducing 
modifications to this kind of models in such a way that 



more realistic sequences of events are obtained. The two 
crucial modificationshav we have proposed are^^ (see a 
more detailed account in the next section) the considera- 
tion of variable thresholds (instead of the constant value 
considered by OFC), and the incorporation of a mecha- 
nism of internal (also called structural) relaxation, quan- 
tified by a parameter R. With these two ingredients we 
have obtained a model in which: 1) the spatial and tem- 
poral correlation of events is comparable to real ones, in 
particular, aftershock sequences obeying the Omori law! 
are obtained, 2) the avalanche size distribution has a GR 
form, with an exponent b very close to actually observed 
values (around 1), 3) the friction properties derived from 
the model reproduce non-trivial results such as velocity 
weakening or the stress peak in experiments of slip-stop- 
slip4 It is worth mentioning that as long as R is larger 
than some threshold value, all these results are obtained 
irrespective of the precise value of R. 

The number of realistic features that have been ob- 
tained with this model, encouraged us to go further and 
investigate in more depth its properties to better under- 
stand how they originate in the existence of the relax- 
ation process. This is the aim of the present work. To 
make the presentation self contained, in the next section 
the model is explained with all necessary detail, in the 
case without relaxation R = 0, and in the presence of 
relaxation R ^ 0. Also, a limiting case of infinite relax- 
ation R — > oo is introduced, that will be interesting to 
consider along the paper. 

In section III we review some results obtained in the 
case R = 0. It is shown that this case is equivalent 
to the problem of depinning of an elastic interface in a 
random medium. We give also a comparison of the re- 
sults obtained at "constant driving velocity" with those 
at "constant force" , showing that there is a clear mapping 
between the two. 
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In Section IV we show that when R ^ we cannot 
make an equivalence between simulations at constant ve- 
locity and constant force. In fact, in constant velocity 
simulations at R ^ 0, there are systematic fluctuations 
of the stress on the system that make the problem non- 
stationary, and this generates a scenario in which a de- 
scription in terms of a "seismic cycle" naturally emerges. 

In Section V we discuss the scaling of the linear size 
of the largest events L max with the control parameter a. 
We show that with an appropriate interpretation of a 
as inversely proportional to the thickness of the sliding 
slabs, there are events that are comparable in size with 
the size of the system. This is a consequence of the ex- 
istence of relaxation, and does not occur in the R = 
case. 

Section VI is devoted to the following issue: the model 
with R = gives a consistent (although unrealistic) value 
for the b exponent of about b ~ 0.4. The inclusion of 
relaxation drives this exponent to values closer to the 
realistic value b ~ 1, without tuning of any parameters. 
What is the reason of this behavior? As we will see, we 
are only in the position to provide a partial answer to 
this question. 

Finally, in Section VII, we make a brief summary and 
present our conclusions. 



II. MODELS 
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A. OFC model 



FIG. 1. A sketch of the sliding situation that we are studying, 
(a) Two solid blocks slide against each other due to a constant 
driving between the top and bottom planes. Their relative 
velocity is V. Dimensions of the blocks are L x , L y in the 
sliding plane, and L z perpendicularly to it. (b) The solids 
in (a) are replaced by a rigid surface and an array of small 
blocks joined by springs with stiffness ko- Driving acts on 
each block through a spring of stiffness k\. Note that the 
value of fei is inversely proportional to L z in (a), (c) The 
values of friction forces Ui and maximum forces u\ at each 
site are shown for the OFC* model. Friction forces increase 
uniformly (as indicated by the arrows) as driving proceeds, 
(d) The prescription as one m reaches the local threshold: 
Ui decreases by one unit, and the neighbors increase their u 
values a quantity a. 



Our starting point is the cellular automata model pro- 
posed by Olami, Feder, and Christensen^ (OFC) to de- 
scribe seismicity. The OFC model considers a set of real 
valued variables Uj where i indicates the position in a 
two dimensional lattice. Ui is interpreted as the force 
that a rigid substrate exerts on a solid block at position 
i, and it represents the local stress between the sliding 
plates (see Fig. [1}. The system is driven by uniformly 
increasing the values of Ui with time at a rate V, simu- 
lating the tectonic loading of the plates. Every time one 
of the variables Uj reaches a maximum value (ordinarily 
set to an uniform, dimcnsionlcss value of 1), the local 
stress Ui is 'discharged' by setting it to zero. This local 
stress drop Au produces a stress increase onto neighbor 
blocks according to Uj — > Uj + aAu, where j indicates 
a neighbor site to i. The value of a can vary between 
and a c = 1/Z, Z being the number of neighbors in the 
lattice. We will refer only to the case of a square lat- 
tice, so Z = 4, a c = 1/4. The case a = a c is called the 
conservative case, whereas a < a c are non-conservative 
cases. In the mechanical interpretation of the model, the 
value of a is related with the stiffness of the springs that 
inter-connect the blocks, and springs that push from the 
blocks at velocity V (see Fig. [TJb)). The actual relation 
is a = k /(ki + 4fc ). A discharge event can produce the 
overpass of the maximum local stress on one, or more 
than one neighbor and in this way a large cascade can 



be generated. This cascade is called an event, and is 
identified with an individual earthquake (note that the 
complete cascade is assumed to occur at constant time, 
namely, earthquakes are instantaneous). The size S of an 
event is defined as the sum of all discharges that compose 
the event, and the magnitude is defined as M = | log 10 S, 
so to match (up to an additive constant) the usual defi- 
nition used in geophysics^. 

The OFC model is typically simulated using open 
boundary conditions, otherwise the spatial homogene- 
ity generated by the use of periodic conditions induces 
a strong global synchronization in the model. The OFC 
model displays an exponential decay of number of events 
as a function of magnitude (or a power law decay if ex- 
pressed in terms of S) compatible with a GR law of the 
form N(M) ~ 10~ bN . The 6-value depends on the value 
of a. Realistic values of b are obtained for a ~ 0.2. 



B. OFC* model 

The maximum values that the variables Uj can with- 
stand in the OFC model, are set to a constant, uniform 
value of 1. Having in mind a realistic situation of a het- 
erogeneous fault, with the constitutive materials having 
different properties at different positions, it becomes nat- 
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ural to consider a case in which the threshold values are 
not constant but have some spatial variation. This is the 
first modification we introduce into the OFC model. In 
concrete, the values of the local thresholds are called uf 1 , 
and we draw them from some random distribution. Each 
time Ui overpasses the local threshold u 1 / 1 , m is updated 
to a new value. We use the update rule Ui — > m — 1, 
i.e, we implement a prescription of a unitary local stress 
drop. Upon this drop of the local stress, the values of u on 
neighbor sites are updated as before, namely Uj — > Uj+a, 
for j neighbor to i (see Fig. [TJc-d)). Every time Ui is up- 
dated, a new value is assigned to the local threshold u\ , 
taken from its original random distribution. This pre- 
scription is justified on the same physical arguments as 
before, since the sliding pieces can reasonably be thought 
as finding different maximum strengths as sliding pro- 
ceeds. We refer to the OFC model with this modification 
as the OFC* model. 

The results obtained with the OFC* model differ qual- 
itatively from those obtained with the OFC model. In 
particular, the b exponent becomes a- independent, tak- 
ing a value around b ~ 0.4 (yet quite unrealistic for earth- 
quakes). Interestingly, now the maximum avalanche size 
is controlled by the value of a, diverging for a — > a c (as- 
suming an infinite spatial extent of the system). This 
may be considered a weak point, since it seems to in- 
dicate that is not sufficient to have an infinitely large 
system to have critical behavior, the condition a — > a c is 
also needed. However, we will see in Section V that this 
argument is flawed, as the large system size limit also 
implies L 2 — > oo, and in this limit we naturally obtain 
that a — > a c . 

C. OFCR model 

The second and crucial modification we made on the 
OFC model is the introduction of internal relaxation ef- 
fects. We will not repeat here the arguments that justify 
its introduction, but let us just mention the following fact 
that points to the necessity of these effects: In the OFC, 
or OFC* models, a sudden stop of the tectonic loading 
makes all seismic activity to cease at once. This is un- 
realistic, as for instance the processes that trigger after- 
shocks after a main event depend on the rearrangements 
that occur in the fault after the main shock, and are 
not directly related to tectonic loading. In other words, 
dynamical effects within the faults should be taken into 
account in order to have a realistic description of the 
seismic process. This is what we did by introducing the 
internal relaxation mechanism into the OFC* model. As 
this relaxation is controlled by a parameter R, we will 
call this the OFCR model. The concrete implementation 
prescribes that the evolution of the variables Ui is given 
by 

^ = R(V\) i + V (1) 



namely, in addition to the external loading (represented 
by the V term) there is a tendency to make the values 
of Ui progressively more uniform in the system. This 
relaxation term produces the appearance of aftershocks 
and realistic friction properties in the system. In addi- 
tion, the exponent of the GR law is modified, acquiring 
a value consistent with observations.^ 

For the rest of the paper it will be important to con- 
sider a couple of variations of the previous relaxation 
mechanism. One of them is to consider relaxation in a 
"mean field" manner. By this we understand that in- 
stead of the sum over the neighbor sites implied by the 
Laplacian operator in Eq. (fTJ) we take the average over 
the whole system, i.e, noting u the mean value of u,;, we 
have 

^ = -4ft(m — u) + V (2) 
at 

(the factor of 4 is included to have a more direct compari- 
son between the R parameters in the two equations). The 
consideration of relaxation in mean field is important be- 
cause it allows a very efficient determination of the next 
site that becomes unstable (i.e. the first site for which 
Ui reaches uf 1 ), as Eq. @ can be solved analytically, 
whereas the same determination for Eq. (TTJ) requires a 
time consuming integration scheme. Results using both 
schemes of relaxation do not differ substantially. A com- 
parison between the two schemes was presented in [5] . 

The second special case to be considered is the limit of 
"infinite relaxation" defined by the condition R/V —> oo. 
In this limit, there is a clear cut separation between 
events that are triggered by the tectonic loading V term 
in Eqs. (JXJ) or and those that are triggered by the 
relaxation term proportional to R. In fact, in this case 
events separate naturally into "clusters" . In between suc- 
cessive clusters, the spatial distribution of it, is uniform, 
and the V term acts, moving up all Ui uniformly until the 
lowest u th is reached. At this point, an event is triggered 
that leaves the it's in a non- uniform state, over which 
the relaxation term acts. Under this action, aftershocks 
are triggered until the spatial distribution of u's becomes 
uniform again, and the process is repeated. 

III. BACKGROUND MATERIAL ON THE R = 
CASE 

We present in this section a few reference results for 
the case of no relaxation, but with variable thresholds, 
namely, what we have called the OFC* model. 

A. Mapping between elastic interfaces in random 
media and the OFC* model 

First of all we show that the OFC* model can be 
mapped onto the problem of a two-dimensional clastic 
interface close to the depinning transition. 
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Let us consider an ensemble of block positions hi iden- 
tifying a two dimensional interface. The energy of a given 
interface can be written as the sum of three contributions: 

(a) the clastic energy between nearest neighbors 
blocks: Ei = ^J2{i fiQ 1 * ~ ^j) 2 , where h i: hj are the 
positions of the blocks i and j, and kg is the "surface 
tension" of the elastic interface. 

(b) the energy of tectonic loading: E 2 = J2i iH^ — 
Vt) 2 , where V measures the rate and ki is the stiffness 
of the loading. 

(c) the random energy associated with the block posi- 
tion hi, that is responsible for the pinning of the blocks. 

Elastic interactions and tectonic loading generate a 
force Fi = —dhiiEl + E2) which competes with the pin- 
ning force F p i nn (hi) induced by the random energy The 
force Fi acting on each block increases with time because 
of the loading. When Fi overcomes F P i nn (hi) , the block 
i becomes unstable and its position is increased by one 
(hi — > hi + 1). As a consequence (i) the force Fj drops to 
Fi — (4fco + ki), (ii) a new pinning force is associated to 
the new position hi + 1 and (iii) the force acting on the 
neighboring blocks increases Fj —> Fj + k Q . 

The dynamics of this interface in presence of a disor- 
dered landscape can be easily identified with OFC* dy- 
namics in presence of variable thresholds, the stress Ui 
being Fi/(ki + 4/co) and the control parameter a being 
fc /(fci +4fc )- 



B. Review of the depinning transition 

The zero temperature motion of an elastic interface 
in random media has been studied in detail in the last 
two decades^—. Two different dynamical protocols are 
possible: (i) constant velocity: which corresponds to 
the OFC* model with a tectonic loading, or (ii) con- 
stant force: where the parabolic tectonic loading, E%, 
is replaced by the action of a constant and uniform 
stress a, and the force acting on the block i writes 
Fi = {-d hi E x )+Ak a. 



1. Constant force 

For a given (and low) a the interface is pinned in a 
metastable state for which the forces acting on all the 
blocks arc smaller then the forces induced by the pinning 
centers. When a is increased of an infinitesimal amount 
5a two things can happen: (i) either the metastable state 
remains metastable, or (ii) the metastable state becomes 
unstable and moves to a new metastable position. The 
distance between the old pinned position and the new 
one can be measured as the volume S included between 
the two consecutive metastable states. This volume rep- 
resents the size, S, of the event. A critical threshold, a c 
exists, above which there are no metastable states, and 
the interface evolves indefinitely. The transition between 



the pinned phase, a < a c , and the moving phase, a > o~ c , 
is called depinning. 

Results of simulations at different values of a < a c 
show that the distribution of events follow a power law, 
with a large size cutoff, S ma _ x (a) and an exponent r > f , 

N(S)<x±f(—?—^. (3) 

The exponent r used by the statistical mechanics com- 
munity is related to the exponent b typically used by 
geophysicists when looking at real earthquakes. The re- 
lation between them is given by r = f + 2b/3. We will 
refer mainly to values of r, from now on. The function 
f(S/ S max ) has a fast decay when S ^> S max and the char- 
acteristic large size, S max (a) diverges as a reaches a criti- 
cal threshold a c . The divergence of Smax goes with the di- 
vergence of the linear size of the maximal avalanche L max - 
L max is often seen as the "growing correlation length" as- 
sociated to this dynamical phase transition and, close to 
fc! imax behaves as (a c — a)~ v where, the positive expo- 
nent v is defined in analogy with the correlation length 
exponent of equilibrium critical phenomena. 

Large avalanches occur because the system is "orga- 
nized" on large distances, another consequence of this 
"organization" is that the interface is a correlated object 
displaying a power law roughness. A positive roughness 
exponent ( can be defined from the growth of the distance 
between blocks |/ij — hj\ ~ [i — J | • An argument based on 
a symmetry of the system (statistical tilt symmetry)^— 
allows to relate the exponent v to the exponent £, yield- 
ing 



From basic dimensional arguments we can write 

5 max ~ L 2 J< ~ (a c a)~^ 2+ ~ (a c a)~& (5) 

An important scaling relation between r and £ can be 
established if we compute the average size of an event 5*. 
On one side, using the scaling form of Eq. [31 this quan- 
tity can be written in terms of the cutoff size, in the form 
S oc S^J, as long as 1 < t < 2 which is always the case 
here. On the other side we observe that the displace- 
ment of the center of mass xqm of the interface when 
the stress in increased by an infinitesimal amount can be 
written as, x C m(o- + 5a) - x CM {a) = SN a (Sa)/(L x L y ), 
where N a (Sa) is the number of events when the stress 
jumps from a to a + 5a. If this number does not diverge 
when the threshold is approached (as confirmed by all 
simulations), then we can power expand, and to leading 
order we have N a (5a) ~ ci5aL x L y + . . .. Thus, it fol- 
lows that the average size of an event is proportional to 
S ~ 5xcyi/5a = x, the susceptibility of the interface. As 
the critical threshold is approached the interface moves 
more and more, the position of the interface grows as 
xcm(o~) ~ ^Lax an d thus the susceptibility diverges as 
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X ~ { a c — cr ) We thus obtain 

i/(2 + 0(2 -t) = i/C+1 => 



C + 1 
T = 2 'JTt- (6) 



Using the statistical tilt symmetry relation (Eq. 2]) we 
obtain 



r = 2 



2 + C 



(7) 



2. Constant velocity 

At constant velocity, instead of a uniform stress, we 
have the parabolic potential E2 centered at Vt. The in- 
terface, embedded in this potential can be rough only 
at short length scales. The crossover between the short 
distance regime characterized by a rough interface and 
the long distance regime with a flat interface can be de- 
termined by dimensional analysis of the term E\ ~ L 2< ^ 
against the term E2 ~ kiL 2 ^ +2 . The crossover length can 
be identified with the correlation length L max ~ \j\fk\. 
Using the statistical tilt symmetry of elastic systems it 
is possible to show that this identification is actually 
correct. It is easy to verify that a c — a ~ k\, so that 

«c — Oi — i'max- 



C. Results for the OFC* model 

The previous arguments are consistent with the results 
of our simulations in the OFC* model. In Fig. [2]we show 
a sequence of individual events and the corresponding 
stress as a function of time for a fixed a, and the event 
size distribution obtained for different values of a. A 
value t ~ 1.27 is systematically observed. Also, from 
the distributions of N(S) and the corresponding ones 
for N(L) (L is the typical linear size of the avalanches) 
we obtain the value of the large scale cut off S max and 
the maximum linear size L m ax as a function of a for the 
OFC* model. Results are presented in Fig. [3J We found 
Snmx ~ ipc c — a)~ 7 with 7 ~ 1.37, which is consistent 
with the ratio 7 = 1 + £/2 predicted by the theory of 
elastic interfaces of the previous section. We also find 
from Fig. [3]that L max ~ (a c — a)^ 1 / 2 , showing that sta- 
tistical tilt symmetry is satisfied by the OFC* model, as 
the mapping to the elastic interface problem could have 
anticipated. 

For a fixed value of a we can follow the value of the 
average stress in the system (see Fig. [2]), a = 11%. The 
stress settles around a mean value, with fluctuations typ- 
ical of finite size effects. According to arguments in 
the previous section, the time averaged a vs. a is ex- 
pected to have a critical behavior as a — > a c of the form 
(cr c — cr) ~ (a c — a) 1 ^ 2 ^ for some cr c . This is verified in 
our numerical simulations (Fig. @|, as \/{2u) ~ 0.62. 

The existence of a well defined average value of a for 
each value of a allows to recover the scaling relation 
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FIG. 2. Magnitude-time and stress-time plot for the OFC* 
model, and the event size distribution obtained for three dif- 
ferent values of a as indicated. These distributions are very 
well fitted by a power law with an exponential cut off. 
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FIG. 3. Scaling of (a)S ma x and (b)I/ ma x with (a c — a) for the 
OFC* and the OFCR models, shown in full and open symbols 
respectively. 



(Eq. (J7J) between r and £ using a slightly different ar- 
gument. If an event of size S occurs in the system, it 
produces a reduction of a by an amount proportional to 
SAa/(L x L y ) (Aa = a c — a). If it is assumed (as it is in 
fact numerically verified) that the time interval between 
events does not become singular as Aa — > 0, then there 
is a typical stress increment between events of the order 
of l/(L x L y ). To obtain a stationary mean value of cr, 
these two quantities have to compensate, on a temporal 
average, namely, we can write SAa ~ 1 where with S we 
note the average value of the events. Putting this esti- 
mation together with S oc and using also Eq. ([5]), 
we re-obtain Eq. 0. 

Simulations of the OFC* model at a < a c and con- 
stant velocity are associated (in the mechanical analogy 
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FIG. 4. Time averaged values of a as a function of the values 
of a for the OFC* model. The power law dependence close 
to a c — 0.25, tj c — 0.84236 is emphasized in the logarithmic 
plot in the inset. 



FIG. 5. Magnitude-time plot in the OFC* model, taking V = 
0, a — 1/4, and using the stress a as a control parameter, 
and event size distribution obtained for three different values 
of a as indicated. Note the similarity with Fig. [2] 



of the model) to driving by a finite stiffness spring at a 
constant velocity. The univocous relation between a and 
a allows also an alternative interpretation of the model, 
which is equivalent to driving the model at constant force. 
Instead of applying a finite driving at fixed a < a c , we 
set a = a c , the driving velocity V = 0, and fix the mean 
value of a from outside. A random site is chosen to ini- 
tiate the avalanche. After the avalanche is exhausted, 
the process is repeated. Note that a does not change 
in this process, and remains equal to the externally im- 
posed value. In this way of doing the simulation, after a 
transient period, the avalanche size distribution becomes 
stationary in time, and its statistics coincides with that 
obtained fixing a non-critical a and having a finite V, 
as can be observed in Fig. [5] The expected behavior of 
Smax as a function of a is given in Eq. ([5]). This behavior 
is consistent with our data for which v{2 + Q ~ 2.2. 



IV. OFCR MODEL AND THE SEISMIC CYCLE 

In this section we begin to analyze results obtained in 
the presence of relaxation. We will work in the limit of 
infinite relaxation, and use the mean field implementa- 
tion of relaxation. A typical sequence of events for this 
case is shown in Fig. [SJ We also indicate in that figure 
the evolution of a over the simulation. In the infinite 
relaxation case we will refer to an "internal" and an "ex- 
ternal" temporal scale, which are totally decoupled. The 
external time scale is the scale of tectonic loading, and 
internal time scale is the scale of relaxation in Eq. ((T|). 
In Fig. [5] we plot events as a function of the external 
time scale. This means that at a fixed value of the hor- 
izontal coordinate, we find a whole "cluster" of events. 
The events within some of these clusters, depicted as a 
function of the internal time scale are shown in Fig. E£a) . 
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FIG. 6. (a) Magnitude-time and stress-time plot for the OFCR 
model. The probability of large events is strongly enhanced 
every time some threshold value of the stress (roughly indi- 
cated by the horizontal dotted line) is over passed, (b) The 
corresponding event size distribution for the whole sequence 
at different values of a. 



One may consider a cluster as consistent of a main shock 
and all the aftershocks it produces. The clear cut identifi- 
cation of aftershocks as linked to some main event can be 
unambiguously made only in the present case of full sepa- 
ration of time scales. Note however that within a cluster, 
the largest shock is not necessarily (nor even usually) the 
first one. In any case, we can consider a cluster as the 
full sequence of events that is committed to appear once 
the first event has been triggered, in the case in which 
tectonic loading is stopped after the initial event. We 
emphasize that in the R — > oo case, the value of a shown 
in Fig. |5] is not only the average value of u through the 
system, but it is also the actual value of every Ui. 
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FIG. 7. (a)Magnitude-internal time (U n t) for some of the 
main clusters observed in Fig. [6ja). Unt is measured in units 
of the relaxation parameter, i.e, a whole cluster here occurs 
at a single horizontal coordinate in Fig. |6ja). (b) the event 
size distribution of the individual clusters. 



In the results of Fig. [6] we see a systematic pattern 
of steady increase of a, followed by abrupt drops. We 
interpret this cycle of smooth stress rising and sudden 
stress decrease as the manifestation in our model of a 
"seismic cycle"— . The concept of the seismic cycle sug- 
gests that stress on a fault smoothly builds up for a long 
period of time until is suddenly released in a large earth- 
quake, and the process starts over. It leads to the idea 
that large earthquakes in a given region occur in a quasi- 
periodic manner. This idea, however, has been difficult 
to verify, on one side because of the large time intervals 
between large earthquakes, and also because deviations 
from periodicity could be so large as to transform a quasi- 
periodic sequence in something rather unpredictable. We 
notice however, that there have been a number of at- 
tempts (sometimes successful^) to infer the likelihood of 
a next large earthquake based on the record of previous 
large earthquakes in a particular region. 

In the simulations, an alternation of large events and 
periods of rather low, slowly increasing activity is ob- 
served (see Fig. [6]). However, this sequence is not peri- 
odic, instead it looks like a more or less random process. 
In our case, we have the advantage to have access to all 
variables in the system, particularly to the instantaneous 
value of the stress, and thus we see that a large event oc- 
curs with large probability once some typical value of the 
stress (roughly indicated by the horizontal dotted line in 
Fig. [6]) has been over passed. These observations point 
to the fact that there is some degree of predictability of 
the appearance of large earthquakes in the present model. 
However, we will not discuss here this issue, that we plan 
to elaborate in a forthcoming publication. 

The fluctuations of stress during the seismic cycle for 
the OFCR model are more fundamental than those ob- 
served for the OFC* model in the previous section. In- 



FIG. 8. The seismic moment of individual clusters S c vs. 
the initial value of stress at which they were triggered. The 
crossover stress (corresponding to the horizontal dotted line in 
Fig. [6]) is indicated by the vertical dotted line. When stress 
becomes larger than this value there is some probability to 
trigger a very large cluster (those large clusters have been 
highlighted using a larger symbol size). 



sight into this issue is provided by taking the data in Fig. 
|6] and plotting the total seismic moment S c of each clus- 
ter (namely, the sum of the seismic moments of all event 
within a cluster) which is proportional to the stress drops 
observed in the curve of a(t) as a function of the value 
of stress at which the cluster is triggered. The result is 
shown in Fig. [5] We see that when a is small, the clus- 
ters are typically small also and the stress drop caused 
by these relatively small clusters is not able to compen- 
sate the stress increase caused by tectonic loading. As 
a becomes larger than a certain value (indicated by the 
vertical dotted line in Fig. [8]), there is a finite probability 
that a very large cluster is triggered. The size of these 
large clusters diverges as a — > a c . These large events 
generate the abrupt decrease in stress observed in the 
previous Fig. |6j re-initiating the seismic cycle. This tells 
that the fluctuations of stress in the OFCR model are 
a necessary part of the dynamics of the model, which is 
eminently non-stationary, and the existence of a seismic 
cycle is the manifestation of that. 



V. THE MACROSCOPIC LIMIT OF THE OFCR 
MODEL 

There is a clear cut difference in the macroscopic limit 
of the OFC* and OFCR models that shows that the 
OFCR model is a candidate to be a realistic model for 
describing seismicity, whereas the OFC* model is not. 
This is related with the scaling of the largest events ob- 
served for the two models in the macroscopic limit that 
we intend to discuss now. 

This analysis starts from the key observation that the 
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value of a c — a should be considered as proportional to 
L~ x , where L z indicates the thickness of the slabs that 
slide against each other. This is clear by comparing Fig. 
[TJa) and (b): the value of k\ (and thus of a c — a = 
&i/4(4fci + fco)) is proportional to the restoring force the 
system exerts upon a given displacement of the top plane 
of the slab. For a fixed displacement, the restoring force 
is inversely proportional to L z , namely a c — a ~ L^ 1 . 

In the macroscopic limit we should consider that L 
(L = L x = L y , the two lengths in the x-y plane are 
taken equal for simplicity) and L z go to macroscopic val- 
ues, although not necessarily with the same tendency. 
We will consider a case in which L L z , i.e., we main- 
tain a sort of slab geometry. In this situation we want to 
calculate the maximum spatial extent L max of the earth- 
quakes in the x-y plane for OFC* and OFCR models. 
In the first case the statistical tilt symmetry (verified in 
the numerical simulations, see Fig. [3Jb)) implies that 

imax ~ V( a c - a) 1/2 , from which L max ~ L 1 / 2 . This 
means that as L z goes to infinity, the maximum spatial 
extent of the earthquakes becomes comparatively small. 
In other words, at large scales the sliding is smooth in 
the OFC* model, earthquakes do not survive the macro- 
scopic limit. 

The situation is qualitatively different for the OFCR 
model. Here, numerical simulations indicate that the 
maximum spatial extent of the events has a different de- 
pendence with (a c — a). The results presented in Fig. 
[31(b) indicate a dependence L max ~ l/( a c — a), which in 
terms of L z , can be written as L max ~ L z . This is a re- 
markable result. It indicates that no matter how large L z 
is, the maximum size of earthquakes is a sizable fraction 
of it. We can thus say that in this case, in the macro- 
scopic limit earthquakes persist, and their maximum size 
scales as the system size itself. We think this is a key 
point to claim that this model can describe earthquakes 
in the macroscopic limit. 

The finding of sizable earthquakes with the size L z 
can be correlated with other properties of the model. It 
is accepted in the geophysics community, that systems 
in which earthquakes occur must have an effective fric- 
tion law of the velocity weakening type^^, namely, the 
average friction force must be a decreasing function of 
velocity, at least in a velocity range relevant for the pro- 
cess. Actually, we have shown elsewhere^ that the OFCR 
model does display velocity weakening, and in the origin 
of this behavior is again the existence of relaxation. It is 
immediate to show that such velocity weakening systems 
cannot sustain uniform sliding, and a stick-slip motion 
must necessarily occur for a sufficiently soft driving (i.e., 
for a — > a c ). Now if we assume a generic dependence 
of L max on L z of the form L max ~ L", a value u> < 1 
would correspond to an asymptotically smooth sliding, 
and since this cannot occur, we should have u> > 1. On 
the other hand, it is very unlikely that L max scales as a 
power of L z larger than one, because in a geometry of 
a slab with thickness L z , a static perturbation at some 
point does not have an effect beyond some distance of 



the order of L z . This takes us to the heuristic finding 
that the exponent effectively found in the numerical sim- 
ulations, namely ui = 1 is the most natural result to be 
expected. 

A further question that we asked ourselves concerning 
this point is the following. The OFC* model displays a 
behavior with uj = 1/2, whereas our results for the OFCR 
model indicate uj = 1. How is the transition between 
these two exponents as the precise value of R is varied 
between (OFC*) and a large value (OFCR)? Although 
this is an attractive question to explore, some investiga- 
tion of the possibilities to provide an answer based on 
numerical simulations convinced us that we are not able 
to do so at present. We want to stress however our expec- 
tation that an answer to this question could give insight 
into the puzzling problem of creeping faults^, namely 
a piece of a fault in which there is an almost complete 
absence of earthquakes and the sliding is smooth, while 
earthquakes occur in adjacent segments of the same fault. 
In our view this could be related to an abrupt transition 
between ui = 1 for a "normal" segment, to a value uj < 1 
in the creeping fault, and this abrupt transition could be 
driven by a smooth change in the amount of relaxation 
in the system. Clearly this issue is an open line of future 
research. 



VI. WHY b ~ 1? 

One of the results that surprised us more concerning 
the OFCR model is the fact that it displays a realistic 
exponent b ~ 1 (r ~ 1-67), whereas the OFC* has b ~ 0.4 
(t ~ 1.27). This is surprising because the inclusion of 
relaxation was aimed to generate aftershocks (what it 
does), and not to obtain a "correct" GR law. So we 
believe there is a necessity to provide an explanation for 
this behavior. 

In section IIIB we have derived Eq. [7] linking the 
avalanche exponent r with the roughness exponent £ 
of the events of the OFC* model. The steps followed 
to arrive at Eq. [7] can be done also in the case of 
the OFCR model. The only important difference is 
that, as explained in the previous section, the scaling of 
L max ~ (a c — a)"" has a different exponent for the OFCR 
case (ui = 1) compared to the OFC* case (uj = 1/2). The 
result that is obtained for a generic value of ui is 



This modified scaling relation suggests already a larger 
r exponent in the OFCR case. In fact, assuming only a 
non-negative value for £, we get r > 1.5. To determine 
the value of £ we need to characterize the internal struc- 
ture of an avalanche. Four quantities can be used to this 
scope: (i) the already mentioned size S, which counts 
the total number of discharges, (ii) the area A defined 
as the total number of blocks involved in the avalanche, 
(iii) the "duration" T defined as the number of discharge 



steps that are necessary to exhaust the avalanche (note 
that sites that are critical at the same moment are dis- 
charged in parallel, in a single step), and (iv) the lin- 
ear size L of the avalanche. Similarly to the size S, all 
these quantities are power law distributed up to a cut- 
off which depends on a c — a. In the stationary regime, 
the study of the cut-off divergence, as the threshold a c 
is approached, allows to determine the value of the crit- 
ical exponents like £, v, ... A different approach which 
does not rely on the stationary regime, consists in study- 
ing the behavior of A, T and L as a function of the size 
S. For large S we expect a power law behavior with 
(A{S)) ~ S KA , (T(S)) ~ S KT and (L(S)) ~ S* Ki , where 
(A(S)) is the area averaged over all avalanches of size 
S, (T(S)) is the duration averaged over all avalanches 
of size S and (L(S)) is the linear size averaged over all 
avalanches of size S. For the OFC* model, scaling re- 
lations predict in two dimensions: ka = 5^ ~ 0.73, 
kt = 2^ ~ 0.57 an d «l = 2TC ~ 0.36, where we have 
used the value z = 1.56 for the dynamical exponent. Nu- 
merical simulations on the OFC* model agree with this 
predictions. More surprisingly also the avalanches of the 
OFCR model seem to have the same internal structure 
and the same value of the k exponents. As a first conse- 
quence if the value of £ does not change appreciable from 
the value £ ~ 0.75 in the absence of relaxation, through 
Eq. [51 we predict r ~ 1.64, which is consistent with the 
value observed in the simulations (Fig. On the other 
hand the invariance of £ when including relaxation gives 
a clue on the possible mechanism for the change of the 
t exponent. In fact, if the new value of r corresponds 
to the system being controlled by a different fixed point 
in the parameters space (i. e. a new universality class), 
we should expect that all exponents of the OFCR model 
are different from those of the OFC* model. However, 
as explained by Durin and Zapperii^ in the context of 
a magnetic problem, if the internal (k) exponents of the 
avalanches do not change, it is likely that the new r ex- 
ponent originates as an effect of non-stationarity. If there 
is in the model a parameter that samples a wide distri- 
bution of values, the complete event size distribution can 
be generated as an integration of partial distributions re- 
stricted to particular values of the parameter. If, for in- 
stance, this parameter controls the position of the cut off 
of the distribution, such an integration produces a larger 
effective r exponent, but does not affect the exponents 
ka, kt and kl characterizing the internal structure of 
the avalanches. 

A parameter that widely fluctuates in the OFCR model 
is the stress a. As a controls the position of the cutoff 
of the distribution, an integration over all a values could 
be the responsible of an overall power law with a larger 
decay exponent. If this explanation is correct, the event 
size distribution restricted to some small stress interval 
should display a t ~ 1.27, with a cutoff that is stress- 
dependent. However this is not the case. In Fig. |7[b) 
we have shown event size distributions restricted to the 
events in some of the individual clusters observed in Fig. 
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FIG. 9. (a) A typical distribution of forces u and thresholds 
u th in the absence of relaxation. For u th we show in gray the 
distribution from which each value is chosen, and not the val- 
ues themselves. In the presence of relaxation, the configura- 
tion is modified as indicated in (b): sites for which relaxation 
produces an increase of m can trigger aftershocks. If not, 
the distribution of it continues to be exponential. Sites for 
which relaxation produces a decrease of Ui get a gap between 
u and the lowest possible u , and make these sites have a 
lower probability to be destabilized upon receiving load from 
neighbors. 

[5] (and thus characterized by a single value of a at trig- 
gering). We see for each individual cluster an event size 
distribution with the modified exponent t ~ 1.67, indi- 
cating that a is not the parameter over which the effective 
integration takes place. Some other less obvious scheme 
of "integration over a parameter" must be responsible for 
the change in exponent. 

We have not been able to identify clearly in the OFCR 
model the parameter that produces the effective integra- 
tion. However, we have advanced along this line by do- 
ing the following simplified analysis. Let us consider a 
distribution of values of u and u . The thresholds are 
assumed to have an exponential distribution. The typical 
distribution of u and u th in the absence of relaxation is 
sketched in Fig. 0(a). Note that instead of showing the 
actual values of u , we show the distribution from which 
they are chosen. In the presence of relaxation, the values 
of the forces move as indicated in Fig. H0[b) . In this pro- 
cess, any site for which m < a increases its value. As far 
as the actual value of u th is not reached, the probability 
distribution of u th continues to be exponential above Ui. 
If Ui reaches uf 1 an aftershock is triggered. On the other 
hand, sites for which Ui > a evolve decreasing its Ui value 
because of relaxation. 

To understand the effect of relaxation, and having in 
mind the previous scheme, we do the following simula- 
tion. We take uncorrelated values of u, from some distri- 
bution of width 8u and mean u, and take thresholds from 
an exponential distribution above the corresponding u, to 
mimic the configuration in Fig. EJa). Then, the values of 
u are shrunk towards a by a factor s (0 < s < 1) to sim- 
ulate the effect of relaxation, namely u a-\-s(u — a). If 
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FIG. 10. Avalanche size distribution for an uncorrected dis- 
tribution of it's and a shrink parameter s, as indicated. The 
distribution of u's is fiat between u ± Su, with u = 1.29, 
and Su — 1. The result of an integration over the parameter 
< s < 1 is shown by the line labeled "averaged" (vertically 
shifted, for convenience). As a reference, power laws with 
exponents 1.27 and 1.67 are also plotted. 



in this process u becomes larger than the corresponding 
u , the threshold is chosen again. On this final config- 
uration we destabilize a site at random and measure the 
size of the avalanche that is generated. Statistics is col- 
lected on many realizations of the process. Results are 
shown in Fig. [TDJ For s = 1 (i.e., no relaxation) we get a 
power law distribution with an exponent close to (actu- 
ally, a bit larger than) the normal value r ~ 1.27. As s is 
taken lower than one, the distribution becomes steeper. 
The qualitative reason of this lies in the gaps that appear 
between the u values and the corresponding thresholds 
(Fig. H^b)): sites with larger gaps become less probable 
to be activated upon receiving load from a neighbor. 

By choosing s from a flat distribution between and 
1 (i.e., doing an integration over the value of s), a power 
law decay is recovered, with an exponent similar to that 
found with the full OFCR model. In a full simulation 
with the OFCR model, one can imagine a situation in 
which each event occurs at a particular value of s, and 
an effective integration over this parameter occurs, gen- 
erating the modified exponent. A full justification of this 
scenario is not trivial, and we expect to provide it else- 
where. In any case, our conclusion is that an effective 
integration over an internal parameter that is related to 
relaxation is the responsible for the change of the t ex- 
ponent in the OFCR model. 



VII. SUMMARY AND CONCLUSIONS 

In this paper we have elaborated upon the properties of 
a seismicity model that is based on the one proposed by 
Olami, Feder, and Christensen, and that incorporates re- 
laxation effects as a fundamental ingredient. In the past, 
this model was shown to generate realistic sequences of 
events, in particular displaying aftershocks following an 
Omori law. It also generates a Gutenberg-Richter event 
size distribution with realistic power law behavior, and 
reproduces frictional properties experimentally observed 
in geological materials. 

Here we have focused on the following issues. First of 
all, we have shown that there is a natural notion of a 
"seismic cycle" in the model. This occurs because the 
stress in the system has temporal variations that reflect 
the overall state of the plates, and we have found that 
mayor earthquakes can be expected only when this stress 
is larger than some minimum value. Second, we have 
made an analysis of the scaling of the largest events in 
the system upon system size increase. With the appro- 
priate interpretation of the parameter a of the model, it 
was suggested that in fact the size of the largest events 
scales as the system size (in particular, with the thickness 
of the plate we are trying to model) . This indicates that 
the model is able to describe earthquakes in the macro- 
scopic limit. Note that this is not true in the absence of 
relaxation, since in this case the largest events scale with 
system size with a lower-than-one power. 

Finally we have tried to understand why relaxation is 
able to tune the exponent of the GR law to a realistic 
value around 6 = 1 (or r ~ 1.67). We provided evi- 
dence that the reason of this modification is mainly the 
creation of "gaps" between the u values and the corre- 
sponding thresholds u th , originated in the existence of 
relaxation. Despite this qualitative understanding of the 
effect of relaxation, we have not been able to explain why 
the exponent of the GR law in the presence of relaxation 
seems to adjust systematically around the value b ~ 1 
(or r ~ 1.67). Taken as a whole, the results presented in 
this paper reinforce our view that the relaxation mecha- 
nism introduced on some seismicity models provides an 
important tool to study seismic phenomena. 
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